Iranian Journal of Numerical Analysis and Optimization 


Vol. 12, No. 3 (Special Issue), 2022, pp 497-512 
DOI:10.22067 /ijnao.2022.73876.1079 
https: //ijnao.um.ac.ir/ 


Fe 


How to cite this article ro} 
BY NC Research Article 


Estimation of the regression function 


AMS 


Keywords: Legendre wavelet; Operational matrix; Wavelet approximation; 


Legendre wavelets 


M. Hamzehnejad* , M.M. Hosseini and A. Salemi 


Abstract 


We estimate a function f with N independent observations by using Leg- 
endre wavelets operational matrices. The function f is approximated with 
the solution of a special minimization problem. We introduce an explicit 
expression for the penalty term by Legendre wavelets operational matri- 
ces. Also, we obtain a new upper bound on the approximation error of a 
differentiable function f using the partial sums of the Legendre wavelets. 
The validity and ability of these operational matrices are shown by several 
examples of real-world problems with some constraints. An accurate ap- 
proximation of the regression function is obtained by the Legendre wavelets 
estimator. Furthermore, the proposed estimation is compared with a non- 
parametric regression algorithm and the capability of this estimation is 
illustrated. 
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1 Introduction 


Let f : [a,b] —> R with independent observations {(#;,y:),4=1,...,N}. 
Consider the following nonparametric regression model to provide an estimate 
for f: 

yi = f(t) +e, (1) 


where x; € [a,b] and e; have Gaussian noise. It is well known that the 
following optimization problem approximate the regression function f [7, 13]: 


: 1 ~ 2 A : Tr 2 

Ry LO Fea + x fC) (2) 
where S denotes the set of functions f satisfying the constraints and the 
constant A is called smoothing parameter. The first term measures close- 
ness to the data, while the second term penalizes curvature in the function. 
This optimization problem appears in many branches of applied mathematics 
including economics, stochastic processes, statistics, machine learning, and 
control theory, and several studies have been conducted to determine the 
function f [7, 9, 18, 5, 13). 

Using linear combinations of basis functions, such as orthogonal polyno- 
mials, wavelets, and splines is a popular approach to estimating the function 
f [7, 18, 5, 17, 11, 3, 6, 16]. This kind of method can be expressed as a matrix 
equation that contains a penalty term. Although it is not possible to get a 
clear and accurate answer to this problem, it is necessary to use approximate 
methods to solve it. Calculating the penalty term is an important issue for 
the authors. Wand and Ormerod [18] obtained an exact explicit expression 
for each entry of the penalty matrix by solving numerical integrals. 

It is well known that a single method cannot work for all functions without 
any restrictions. Some of these restrictions include monotonicity, convexity, 
unimodality, or combinations of several types of constraints. For example, 
Mammen et al. [8] considered the regression function under the monotonicity 
constraint and Meyer [9] considered the regression function under constraints 
of convexity and monotone. Also in [1, 12], the authors considered the re- 
gression function under combinations of several types of restrictions. 

In this paper, by using properties of the Legendre wavelets, we provide an 
exact explicit expression for the penalty term only by matrix multiplications, 
which reduce the complexity of the problem. Also, an accurate approximation 
of differentiable functions is obtained by Legendre wavelets. For this purpose, 
we provide an upper bound for the first term of (2). Moreover, by using the 
examples that have been mentioned in [9, 1, 4], we show that the Legendre 
wavelets are a good candidate for the estimation of regression functions under 
various constraints. 

The rest of this paper is organized as follows. In Section 2, we state some 
definitions and properties of the Legendre wavelets. Furthermore, we recall 
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the operational matrix of derivatives, and by using this operational matrix, 
we provide an exact explicit expression for the penalty matrix. In Section 
2, a new upper bound on the approximation error of the partial sums of 
the Legendre wavelets is presented. In Section 3, the performance of the 
proposed estimation is compared with a nonparametric regression method, 
by numerical examples. 


2 Legendre polynomials and wavelets 


In this section, we study Legendre polynomials and wavelets by presenting 
some necessary definitions and theorems. The well-known Legendre poly- 
nomials are defined on the interval [—1,1] and can be determined by the 
following recurrence formulas [15]. 


(m + 1) Lm-+1(#) = (2m a 1)¢Lm (2x) _ mLm—1(2), m= 1, 2, 3, avai) 


where Lo(t) = 1 and Li(a) = x. The following relation is hold for Legendre 
polynomials [15, eq. 3.176a] 


Lal) = = (Eins (t) ~ Lona) (3) 


Moreover, we have the following uniform bound for Legendre polynomials 
[15] 
\Em(a)| <1, « €[-1,1], m>0. (4) 


Legendre wavelets are defined on the interval [0, 1] as follows: 


Un m(t) = (m+ iy or (yd), Berea 
| otherwise, 


where k € N, m=0,1,...,M—1, and n = 0,1,...,2* — 1. The Legendre 
wavelets are an orthonormal basis for L? (0, 1] and the following orthogonality 
holds: 


1 
| Vm,n (t)pr,5(t)dt = dmr Ons: 
Let f(t) € L? [0,1]. Then 


2k_1 M-1 


f(t) =~ S- » Cnym¥n,m (E) = CTW(t), 


n=0 m=0 


where Cn.m = i f(t)tnm(t)dt. The vectors C and W(t) are 2*M x 1 vectors 
given by 
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C = [co,0, see »CO,M-1; C1,05 pare C1,M-1), sey C2k-1,05 sees Cok_1,M—1| 


5] 


W(t) =[vo0(0),---, Yo,m—1(t), Y1,0(t),---, Y1,m—1(t), 
-+,Por_r0(t),--- wee) 


The Legendre wavelets approximation finds a shape constrained f to the 
minimization problem (2). In the minimization problem (2), we set 


2k_1 M-1 


D2 > Gala) 


n=0 m=0 


For simplicity, we can set¢(j—1)x m)+j41(t) = Wig (t) and e(@—1)xm)+j41 °= 


ey fort = 1,..., 2° and j = 0,...,M— 1. Hence the following vectors are 
obtained: 
Pp T 
W(t) = [di(t),..-, Poem (t)] ’ C = |c1,¢2,...,Coem] . (5) 

Therefore, we have 

2*M 

FE) = Do evs (0), 
j=l 


where w;(t) are the Legendre wavelets. Therefore the objective function to 
minimize (2) is the following penalized least square: 


1 2*M 3 M ‘ 
min WN » Yi — > C5 tb; (vi) Tae? Hee dt, 
w=1 j=l 

where 

1 [2*M 2* M 2k M 

[ [Sowa] a= T Yas [ weowrna 
0 \ j=1 i=1 j=l 
N 

Suppose that V is a matrix by elements of the form Vij := + » Wi (ar); (21), 


that P is a matrix by elements Pj; = ie ys”) (y(t )dt, Ae that the ele- 
N 


ments of vector b are defined by b; = x > wi(ai)m, i,j =1,...,2*-1M, so 
I=1 

the minimization problem (2) has the following quadratic form of minimiza- 

tion [5]: 
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1 1 

min =O? VC—bC + (=C7PC). (6) 
CeR?*M 2 

By taking the derivative of (6) in terms of C' and put it equal zero, we obtain 


the following equation: 
(V+ AP)C =b. (7) 


Now focus on the second term, to determine an appropriate operator matrix 
to solve the problem (2). An important issue is to calculate the elements of 
the matrix P. We use Legendre wavelets operational matrix of derivative, 
to get the new structure of the matrix P. The following theorems determine 
the Legendre wavelet operational matrices of derivatives, which are used to 
solve differential equations. 


Theorem 1. [10, Theorem 1] Let W(t) be the Legendre wavelets vector as 
in (5). Then the derivative of the vector Y(t) can be expressed by 


dW (t) 
— = Dut), 


where D is the 2*M operational matrix 


Ce ree 
D= 


where F is an M x M matrix such that (r,s)th entry of F' is defined as 
follows: 


r=2,...,M, 


dd 
s=1,...,r—1, ecpe/oed, 


ee grrr —1)(2s — 1), 


0, otherwise. 
Theorem 2. [10, Theorem 2] By using Theorem 1, the operational matrix 
for nth derivative can be derived as 


d"U(t) 
di” 


= D" W(t), 


where D” is the nth power of the matrix D. 


Therefore, using these operational matrices, the elements of the matrix P 
in (7) are introduced in the next theorem. 


Theorem 3. Let U(t) be the Legendre wavelets vector defined in (5). As- 
sume that r is a nonnegative integer and that the elements of the matrix 
P= [P| are 2 = i. wh” (tO (e)de. Then P;; has the following exact 
explicit expression 
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-_ r r\T ee k 
Pj = (D;)(D5) ’ i,j =1,...,2 M, (8) 
where Df is the 7th row of the operational matrix D” as in Theorem 2. 


Proof. By using Theorem 2, the elements of the matrix P are as follows: 


Pi; -[ W(t wi(na= | wrw@y(Dpwayat i,j =1,...,2°M. 
0 


(9) 
2M 
Let D?w = d\.(t). Then 
Pi -[ (Pui) + aa +d) Yoru (t)) (dvr) + $d. Horas (t)) dt 
,2*u aku 2k Moku 
=| > ye a ws (tyr(eat = S> >> dp fv = (t) Wy (t)dt 
s=1 [=1 s=1 [=1 
According to the property of orthogonality, we have 
1 
| wWs(t)uy(t)dt = dgt. (10) 
0 


By using (10), -> ai) i= Di (Diy. 


Therefore, we can calculate the elements of the matrix P only by a matrix 
multiplication. By solving system (7), the appropriate weight coefficients are 
obtained to approximate the function f. 


3 Error analysis 


In this section, we present an error estimate of the partial sums of Legendre 
wavelets to the regression function f. For this purpose, we benefit from the 
well-known mean-square error (MSE). By using the MSE [16], we measure 
the performance of the estimator f as follows: 


2 
MSE(f,f) n> E fic f(xi)| 
The Legendre wavelets estimator f can be written as 
: ; 2*_1M-1 2k_1 M-1 
f =(f(21),---, 2 D, aebles: rp De creatine) . 
n=0 m= n=0 m= 
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We present a new approximation error of the function f, using the partial 
sums of Legendre wavelets. We know that 


f= S- > CnymPn,m(t) 


n=0 m=0 
2*_1M-1 2*_1 00 oO 00 

= S- Si Cana) - S- S- Ca ate) + S- oy Chi Vinal ts 
n=0 m=0 n=0 m=M n=2k m=0 


(11) 


The last part in (11), S- a CnmWn,m(t) = 0, because the Legendre 
n=2k m=0 
wavelets np m(t) are zero outside of the interval [0,1]. Then 


2*-1M-1 ‘ 2°-1 oc : 
f(t) ad S- ie Cn,mWn,m (Et) = S- S- Cn,mWn,m (Et) 
n=0 m=0 n=0 m=M 
2°-1 oc 
<2 DY lena? llebmm (@)I0- 
n=0 m=M 


We know that ||wnm/(t)||? = 1. Therefore 


: 9 2k_1 M-1 - 2F_1 
F)- FO] =F - XY combamOl] < XY lewml?- 02) 
n=0 m=0 n=0 m=M 


Hence, the approximation error of the truncated series of Legendre wavelets 
depends on the Legendre wavelets coefficients c,,,,. Now, we obtain an upper 
bound for Legendre wavelets coefficients. 


Theorem 4. Suppose that k € N and that f, f’,...,f are absolutely 
continuous on [0,1]. Suppose that V = max {V,, n =0,...,2* — 1}, where 


n+1 


“ok 
Vn= f [Fo*D (| ae, WEG OF a: 


Then form >r-+1, 


M4 
2rk (Qm—2r+3)---(2m—1)(2m+43)++-(2M+2r—1) y/2* (2m—2r+1)’ mode, 
len,m| S V 
,» r even. 
27k (2Qm—2r+3)---(2m+1)(2M+5)--:(2Mm+2r—1)\/2* (2n—2r+1) 
(13) 


Proof. For each 0 <i <r, define the following sequence 
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“2k 
Om =f FO WnamlEat 
1 ae 
=yf(mt 52 f° fO(Lm(2**t—(2n+1))dt, (14) 


where oO, =Cn,m- Let x = 2**1¢— (2n +1). Then 


1. cir ft roi) {@+2n+1 dx 
=\f(m+ 5)2% fs as Lm() 44 
(m+ 3) p e+2n+1 
= ra Ma (a) Lm (a)dx. 


(15) 


By using the equation (3), we have 


(m+ 3) I z+2n+1 
of), = acre fr) (a) (44 @) _ Lin —1(2)) dx. 
(16) 


Using integration by parts, we have 


(m+ 5) 


(nr) (ry (et2n41 : : 
ee 2°" (2m +1) fi ( abet ) aa Lm-1t0)| 4 


ae pL 00 (Re) oneal) ~ ale) 
(17) 


Using the properties L,,(1) = 1” and L,,(—1) = (—1)™ for m > 0, easy 
computations shows that the first term of (17) vanishes. Thus we have 


(r) \/ (m+ 4 fern g+2n+1 
nm = SEF ORT Om + 1) i al 
™m 
(18) 


From (14) and (18), we obtain the following relation between the coefficients 
1 
(Bo we ED 2 HT) 
Cnym — 2k+1(2m a 1) (oy ony) : (19) 


Now, we obtain an upper bound for ve. We can see that 
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fe 1) oa 
Roce) sae A i porn (< +2n+1 
nym k+1 


bat Lm (x)dx 


n+1 
= 4/(m+= pee DOL Ot ne 1)ds: 


From (9) and by easy computation, we obtain 


n+1 
i “ok 
lea = i[tm+ sane f * ors (| |m(2**4t — (2n + 1))| at 
ak 


< (am-+ yak J a |fo+D ( t)) at < V,/2k(2m+1). (20) 


Applying (20) in (19), we have 


i (r+1) | (r+1) ) 
a, 
= 2k+1(2m +1) (|e Cnym— 1 + (c Cn ym+1 


k = k 
Ne (2m —1)+V/2 (2m +3). (21) 


FS +1) 


Ce 


Since 


J2m —14+ V2m4+3 < 2/2m +1, 


(21) becomes to 


k 
on), Eee (2m+1) _ V . (22) 
2+D (2m +1) ~ \/2k(2m +1) 


Also, by using (22) in (19), we obtain the following upper bound for ee a 


or 7 ail) 


x 


ce 1 | ake 


ol » i 
num | = Dk+T(2m + 1) 


IA 


1 V i V 
2k+1(2m +1) \ /2k(Q2m—1)  /2*(2m +3) 


= V (2m + 3) + \/(2m — 1) 
~ ORAL(Qm +1)V2F \ \/Qm— 1)(Q2m +3) 
2 2V /(2m + 3) 
~ Dk+1 (2m + 1)./2*(2m — 1)Qm +3) 
V 


2k (2m + 1)\/2*(2m — 1) 
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If we continue the above process, then by easy computation for an integer 


8s > 2, we obtain the following upper bound for of aa as 


v 

|e) < J) 26-D#(2m—2845)-+\(2m—1)(2mF3)--(2m+2s—3)V/2* (2m—2s+3)’ godd, 

ee a Vv s even. 
2(8-1)k (2m—2s+5)---(2m+1)(2m+45)---(2m4+2s—3)./2* (2m—2s+3)’ 


Then (13) holds when s+ 1 =r. 


Now, we are ready to provide an approximation error of the partial sums 
of Legendre wavelets. We show that if the regression function f is smooth, 
then the partial sums of Legendre wavelets converge to it rapidly. 


Theorem 5. Suppose that k € N and that f, f’,...,f are absolutely 
continuous on [0,1]. Moreover, suppose that Ex, (f(t)) = ie _ fo}. 
Then for M>r+tlandr>l, 


ua r odd 
2(r-1k (2M —2r+41)---(2M—1)(2M+43)---(2M+2r—7)x/2* (2M —2r+4]1)’ : 
Ex,u (f(t) < 4 Vv 
> rT even. 
r2(r—-Dk (2M—2r41)---(2M+1)(2M+4+5)-»:(2M+2r—7)\/2* (2M —2r41) 


Proof. Let r be an odd integer. Applying (13) in (12), we obtain 
Ex,m (f(t) 


Vv 
sd px, 2rk (2m — 2r + 3)--- (2m — 1)(2m +4 3)--- (2m + 2 — 1)\/24(2m — 2r + 1) 


gk 
V ae 1 
ark, /2k(2M — 2r +1 >> D7 (2m — 2r + 3)--- (2m — 1)(2m + 3)--- (2m +4 2r — 1) 
V ee 1 
ark / 2k (2M — 2r + 5% Pa, m + 2r —1)"-1(1 eee Vee (amt arcay) 
2 V 
© 2Pk\/28(2M — 2r + 1)(1 TEMES) sees (aM par=3y) 
2k_4 1 
ij ————_—_—__— dx 
iar ey ae (2a + 2r—1)r-1 
Qkyv 
r2rk(2M — 2r 4+1)---(2M —1)(2M + 3)--- (2M + 2r — 7),/2%(2M — 2r +1) 
Vv 
= 23) 


r2(r—Dk (2M — 2r +1)---(2M —1)(2M 4+ 3)--- (2M 4 2r — 7),/2k(2M — 2r + 1) 


By a similar approach, the results hold for an even integer r and complete 
the proof. 


Remark 1. The aim of this remark is to draw an approximation error for 
a function f(a), using the partial sums of the Legendre wavelets. Consider 
two functions f(x) = 1+ «—0.45exp[—5(x — 0.5)?] and f(x) = ¢x?|c|. The 
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function f(z) = 1+ a — 0.45exp[—5(x — 0.5)?] is infinitely differentiable. 
In Table 1, numerical results are shown for this function for some values of 
M,k, and r. The numerical results obtained from this table indicate that 
by increasing M, k, and r, the partial sums of Legendre wavelets converge 
to the function f(x) rapidly. Also, consider the function f(x) = §2?|2x 


Table 1: Approximation errors of the function f(x) = 1+ 2 — 0.45exp[—5(2 — 0.5)? 
evaluated by Theorem 5. 


Mok r Exm(f(v)) | M ok r  Exm(f(2)) 

10 1 3 1.920x10-° || 10 1 5 6977x10° 
15 2 3 5669x10-° |) 15 2 5 1.686x 10-7 
20 3 3 3.373x 10-6 || 20 3 5 8670x 10-10 


[19]. This function and its derivatives are absolutely continuous on [0, 1] and 
f(x) = |a|. Also, f@ (2) = 2H(x) — 1, where H(x) is the Heaviside step 
function, which is of bounded variation and f(4) (a) = 26(a), where 6() is the 
Dirac delta function. In Table 55, the numerical results are listed for some 
values of M,k, and r. Moreover, the logarithm of absolute errors is displayed 


Table 2: Approximation errors of the function f(x) = 4x?|2| evaluated by Theorem 5. 


Mk or Exm(f(x)) || Mok or Exm(f(z)) 
10 1 3 1067x1077 //10 2 3 1.887x 10> 
15 1 3 3.251x10-° |} 15 2 3 5.747 x 1076 
20 1 3 «1.459x10-> || 20 2 3 2.579x 1078 
in Figure 1. 
10° 1=3, k=2, f(x)=(1/6 2x1 sot 1=5, k=3, f(x)=14x-0.45expf-5(x-0.5)7] 
+ approximation by Theorem 2 + approximation by Theorem 2 


Logarithm absolute error 
Logarithm absolute error 


Figure 1: Approximation error of the functions f(x) = £2|a| and f() = 1+a-— 
0.45exp[—5(a — 0.5)?]. 
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4 Numerical results 


In this section, we present some examples to illustrate the validity and ability 
of the Legendre wavelets. For this purpose, we use some real-world test 
functions. Suppose that (2;,y;),i = 1,...,N are N independent data with 
the same distribution such that X;,7 = 1,...,N have normal distribution, 
that is, v7; ~ N(u,0). Let y; = f(ai) +e; and let x;,¢;, f be independent with 
penalization order r = 2. We consider different kinds of regression functions, 
which have different constraints on interval [0, 1]. 


Remark 2. Choosing the suitable smoothing parameter A is also an im- 
portant issue in solving the minimization problem (2). Corlay [5] showed 
that A = eae is a suitable smoothing parameter, where the quantity oz, 


is the standard deviation, which scales proportionally with x;. Hence, in all 
2r—1 


examples, the coefficient of the penalty term a a wee is used. 


Example 1. Consider two real regression functions f(a) = 15(a — 0.25)? 
[9] and fo(x) = 1+2—0.45exp[—5(x — 0.5)?] [4]. Then f(z) is convex over 
[0,1] and f2(zx) is strictly monotone over [0,1]. Penalized Legendre wavelets 
regression of samples are plotted in Figure 2. 


> Tue curve 7 Tie cave it 
+ data aah fxJ=1+x-0.45*exp(-5*(0.5)) * data 
LW approximation M=4,k=1,N=1000 LW approximation 


f=15(x-0.25) 
M=3,k=1,N=1000 


4 0 sa Ny we HD AN @ © 


oO 01 #02 03 O04 O05 O68 OF O08 OF 1 0 01 #O2 O03 O04 O05 O68 OF O08 O89 1 


Figure 2: Approximate solution for the regression functions f1(#) and f2(2) in Example 
1 


Example 2. Consider the real regression function f3(x) = 15x? sin(3.7x) + 
2 i= exp[(—3(==*)?)] [7, 1], where o = 0.1 and p = 0.3. This function 
is unimodal (first increasing and then decreasing), concave on [0.55, 1], and 
twice differentiable. We approximate the minimization problem (2) for N = 


1000 samples of (2;, y;). In Figure 3, the numerical results are shown. 


Example 3. Consider the real regression function f4(2) = 10(a — 0.5)? — 
exp[—100(a — 0.25)?][4]. In Figure 4, the numerical results are shown. 


In the following example, we compare our method by a nonparametric Re- 
gression (NR) method. NR methods are very sensitive to parameters such as 
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= True curve bd True curve 
si + data ak + data 
LW _ approximation LW approximation 


aL f0)=15x2sin(3.7x)+20/sqrt(20"pilexp(-(x-0.3)0.4)? 


4 4[  £)=15@)sin((3_7)x)+normpdf(x,0.3,0.1) 
M=4,k=2,N=1000 00 


M=4,k=3,N=10 


Figure 3: Approximate solution for the regression function f3(«). 


15 15 
e True curve 1 ° True curve } 
ab {)=10(x-0.5)?-exp(-100(x-0.25)7) + data aL fG)=10(x-0.5)?-exp(-100(x-0.25)7) + 
M=4,k=2,N=1000 LW approximation M=4,k=3,N=1000 LW approximation 
7 


Figure 4: Approximate solution for the regression function f(a). 


the bandwidth selection, the regression order, and the shape of the smoothing 
kernel. In these methods, the choice of order and especially the bandwidth 
parameter can be a hassle [14]. In the previous example, we observed that 
the Legendre wavelets regression (LWR) method provides a good estimate 
for N samples (x;,y;), which does not depend on any parameter except the 
choice of k and M, where k specifies the level of resolution, 2” sub-intervals 
on [0,1], and M determines the degree of wavelets. Note that the selection 
of these two parameters is easy. 


Example 4. Consider the functions fi(z) = 1+ a — 0.45exp[—5(x2 — 
0.5)?], fo(z) = 152? sin(3.72) 4 Ww exp(—$(#=#)?) and f3(x) = —#3 — 2”. 
In Figure 5, we approximate the minimization problem (2) for N = 250 
samples of (x;,y;) and compare this method by a nonparametric regression 
method, which are shown in Figure 5. 


5 Conclusion 


In this paper, Legendre wavelets were used to approximate the regression 
function. A new operational matrix was introduced to simplify the mini- 
mization problem in (2), which is useful for new research in financial math- 
ematics and numerical analysis. Moreover, a new approximation error of 
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Figure 5: Comparing the Legendre wavelets estimation (black curve) with the nonpara- 
metric regression (blue curve). Due to nonoptimal choices of h, under-fitting occurred 
in the first row and over-fitting occurred in the second row for nonparametric regression 
for the functions mentioned in Example 3. 


a differentiable function f using the partial sums of the Legendre wavelets 
was provided. Numerical experiments were performed for a variety of real 
regression functions (see [9, 1, 4]). The proposed method was executed on 
some popular functions, and the numerical results were compared with the 
nonparametric regression method. Finally, the capability of the proposed 
method was illustrated. 
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